%-------------------------------------------------------------------------------

% This file is part of code_saturne, a general-purpose CFD tool.
%
% Copyright (C) 1998-2025 EDF S.A.
%
% This program is free software; you can redistribute it and/or modify it under
% the terms of the GNU General Public License as published by the Free Software
% Foundation; either version 2 of the License, or (at your option) any later
% version.
%
% This program is distributed in the hope that it will be useful, but WITHOUT
% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
% FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
% details.
%
% You should have received a copy of the GNU General Public License along with
% this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
% Street, Fifth Floor, Boston, MA 02110-1301, USA.

%-------------------------------------------------------------------------------

\programme{cs\_boundary\_conditions}

\hypertarget{condli}{}

\vspace{1cm}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Function}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Boundary conditions are required in at least three principal cases:
\begin{itemize}
\item calculation of the convection terms (first order derivative in space) at
the boundary: the calculation uses a flux at the boundary and requires the
value of the convected variable at the boundary when the latter is entering
the domain in the sens of the characteristic curves of the system (in the sens
of the velocity, in the case of the single equation of a simple scalar:
 sufficient interpretation in the current framework of
\CS\footnote{except with the compressible module, see. \fort{cs\_cf\_boundary\_conditions}})~;
\item calculation of the diffusion terms (second order derivative
in space):
a method to determine the value of the first order spatial derivatives
at the boundary is then required
 (more exactly, the terms that depend upon it are required,
 such as the stresses of the thermal fluxes at the wall):
\item calculation of the cell  gradients: the variable at the boundary faces
 are required (more generally, the discrete terms of the equations which depend
upon the gradient inside boundary cells are required, such as the transpose
gradient terms in the Navier-Stokes equations).
\end{itemize}
These considerations only concern the computational variables
(velocity, pressure, Reynolds tensor,
scalars solution of a convection-diffusion equation). For these variables
\footnote{
The other variables
(the physical properties for instance) have a different treatment which will
not be detailed here (for instance, for the density, the user defines
directly the values at the boundary. This information is then stored~; one
is referred to \fort{cs\_user\_physical\_properties} or \fort{cs\_physical\_properties\_update} for more information).
},
the user has to define the boundary conditions at every boundary face
(\fort{cs\_user\_boundary\_conditions}).


The \fort{cs\_boundary\_conditions} subroutine transforms the data provided by the user
(in \fort{cs\_user\_boundary\_conditions}) into an internal format of representation of the boundary
 conditions. Verifications of the completeness and coherence are also
performed (in \fort{cs\_boundary\_conditions\_check}). More particularly, the smooth-wall boundary conditions
 (\fort{cs\_boundary\_conditions\_set\_coeffs\_turb}), the rough-wall boundary conditions (\fort{clptrg})
and the symmetry boundary conditions for the velocities and
the Reynolds stress tensor (\fort{cs\_boundary\_conditions\_set\_coeffs\_symmetry}) are treated in dedicated subroutines.

The \fort{cs\_boundary\_conditions} subroutine
provides as an output pairs of coefficients
$A_b$ and $B_b$
for each variable~$f$ and for each boudary face.
These are used for the calculation of the discrete terms in the equations
to be solved. More specifically, these coefficients are used to calculate
a value at the boundary faces $f_{b,int}$ (localised at the centre of the
boundary face, barycentre of its vertices) by the relation
 $f_{b,int} = A_b+B_b\,f_{I'}$, where $f_{I'}$ is the value of the variable
at point $I'$. $I'$ is the projection onto the centre of the cell
adjoin to the boundary on the line normal to the boundary face and passing
by its centre
(see figure~\ref{Base_Condli_fig_flux_condli}).

See the \doxygenfile{condli_8f90.html}{programmers reference of the dedicated subroutine} for further details.

\begin{figure}[h]
\centerline{\includegraphics[height=8cm]{fluxbord}}
\caption{\label{Base_Condli_fig_flux_condli}Boundary cell.}
\end{figure}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Discretization}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\etape{Notation}
%%%%%%%%%%%%%%%%
In the following, we will denote by {\it VarScalaire} any variable
\begin{itemize}
\item [-] other than the
velocity, pressure, turbulent quantities $k$, $\varepsilon$,
$R_{ij}$, $\varphi$, $\bar{f}$ and $\omega$,

\item [-] solution of a convection-diffusion equation.
\end{itemize}
The denomination {\it VarScalaire} may in particular designate
the temperature, a passive scalar,
a mass fraction or (unless explicitly stated otherwise) the variance of fluctuations
of another {\it VarScalar}. The inferred state variables (mass
volumetric, viscosity...) will not be denoted by {\it VarScalar}.


\etape{Representation of standard boundary conditions in \fort{cs\_user\_boundary\_conditions} \vspace{0,3cm}}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Standard boundary conditions can be provided by
the user in \fort{cs\_user\_boundary\_conditions}. For this, it is necessary to assign a
type to boundary faces \footnote{The assignment of a type is done
by filling in the table
\var{BC\_TYPE}.}. The default conditions are as follows:

\begin{itemize}
\item [$\bullet$] {\bf CS\_INLET}: corresponds to a Dirichlet condition on
all the variables transported (velocity, turbulent variables,
{\it VarScalaires}...), and
a homogeneous Neumann condition (zero flow) on the pressure.

\item [$\bullet$] {\bf CS\_OUTLET}:
        \begin{itemize}
        \item [-] when the mass flow is effectively directed outwards
of the domain, this choice corresponds to a homogeneous Neumann condition
on all the transported variables and\\
$\displaystyle\frac{\partial^2P}{\partial \vect{n}\partial\vect{\tau}_i}=0$, taken into account as
of Dirichlet for the pressure
($\vect{n}$ and $(\vect{\tau}_i)_{i \in \{1,2\}}$
denote respectively the normal vector of the exit face considered
and two normed vectors, orthogonal to each other and in the plane of the exit face).
This condition is explicit using the pressure field and its gradient
at the previous time step.
Moreover, the pressure being defined to a constant, it is readjusted to a point of
output to store the value \var{P0} (this avoids any drift towards very large values
relatively to the maximum pressure difference on the domain)\footnote{When there is not output,
the spectrum of eigen values of the matrix is shifted by a constant value
in order to make the system invertible: see \fort{cs\_matrix\_wrapper\_scalar}.}.
        \item [-] when the mass flow is directed inwards, an
undesirable situation {\it a priori}\footnote{A message indicates
to the user how many exit faces see a mass flow entering into the domain (i.e. a backflow).
},

A Homogeneous Dirichlet condition is imposed on the on velocity (not on the mass flux),
its value downstream of the domain being unknown. The pressure is
treated as in the previous case where the mass flux is
directed outwards of the domain. For variables other than
velocity and pressure, two cases can occur:
        \begin{itemize}
        \item[-] one can impose a condition of Dirichlet to represent the value of the scalar introduced into the field by the faces of boundary concerned.
        \item[-] one can impose, as when the mass flux is outgoing, a homogeneous Neumann condition (this is not a desirable situation, since the information carried on the boundary faces then comes from {\it the downstream} of the local flow). This is the case by default if no value is given for the Dirichlet.
        \end{itemize}
\end{itemize}
\item [$\bullet$] {\bf CS\_SMOOTHWALL}: refer to \fort{cs\_boundary\_conditions\_set\_coeffs\_turb} (or to \fort{clptrg} for rough walls) for a description of the treatment of the boundary conditions of the wall (assumed to be impermeable to
fluid). Briefly, we can say here
that a wall law approach is used to impose the constraint
tangential on speed. The wall may slide\footnote{We must then provide
the components of the wall velocity.}. The {\it VarScalars} are assigned
a homogeneous Neumann condition (zero flow) by default.
To impose a wall value for these variables (for example, in the case of a
wall at imposed temperature) a law of similarity is used
to determine the boundary flux taking into account the boundary layer.
In the case of couplings with \syrthes, \CS
receives a wall temperature and sends a heat flux (decomposed as a pseudo interior temperature
and exchange coefficient to allow relaxation).
The standard pressure condition is a homogeneous Neumann condition.
\item [$\bullet$] {\bf CS\_SYMMETRY}: corresponds to homogeneous Neumann conditions for scalar quantities and to classical symmetry conditions for vectors (velocity) and tensors (Reynolds voltages): see \fort{cs\_boundary\_conditions\_set\_coeffs\_symmetry}.
\end{itemize}

\etape{Representation of specific boundary conditions in \fort{cs\_user\_boundary\_conditions}\vspace{0,3cm}}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
We have seen that the assignment to a boundary face of a standard type
(entrance, exit, wall, symmetry) makes it possible to apply coherent
boundary conditions to all variables in a simple manner for the usual
types of physical boundaries.

It is also possible to define specific boundary conditions in
\fort{cs\_user\_boundary\_conditions}, for each boundary face and each variable
\footnote{Such conditions are direclty defined through the
\var{icodcl} and \var{rcodcl} arrays for each face and variable: examples are provided in
\fort{cs\_user\_boundary\_conditions}.}
(these, like the standard conditions, ultimately amount to mixed-type conditions).

The two approaches are compatible and often combined.
Indeed, the standard (default) boundary conditions
can be overridden by the user where necessary, and the defaults used elsewhere.
In any case, it should be ensured that a boundary condition has been defined for each
boundary face and variable.

Compatibility conditions also exist between the different
variables~(see \fort{cs\_boundary\_conditions\_check}):
\begin{itemize}
\item [-] between, wall, symmetry or free exit, it is important that all the
velocity components have the same type of condition~;
\item [-] when the velocity receives an exit condition, it is important
that the pressure is assigned a Dirichlet-like condition. For more
details, refer to the paragraph relating to the exit condition for
pressure, page \pageref{Base_Condli_Sortie_Pression};
\item [-] when one of the velocity or turbulence variables
is assigned a wall condition, so must all such variables;
\item [-] when one of the $velocity or R_{ij}$ components is assigned a symmetry
condition, it should be the case for all components;
\item [-] when a {\it VarScalar} receives a wall condition, the
velocity must have the same type of condition.
\end{itemize}

\newpage

\etape{Internal representation of boundary conditions}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

{\bf Objective}

The conditions defined by the user are
converted in the form of pairs of coefficients $A_b$ and $B_b$
for each variable~$f$ and boundary face. These coefficients
are used for the calculation of
discrete terms intervening in the solved equation and
make it possible in particular to determine a boundary face value
$f_{b,int}$. It is important to insist on the fact
that this value is, in general, a simple numerical value
which does not necessarily reflect a physical reality (in particular
 to the walls, for the quantities affected by the turbulent boundary layer).
We detail the calculation of $A_b$, $B_b$ and $f_{b,int}$ below.

{\bf Notations}
\begin{itemize}

\item[-] We consider the equation (\ref{Base_Condli_eq_conv_diff_condli})
on the scalar $f$, in which
$\rho$ represents the density, $\vect{u}$ the velocity, $\alpha$ the
conductivity and $S$ additional source terms. $C$ is defined below.
\begin{equation}\label{Base_Condli_eq_conv_diff_condli}
\displaystyle\rho\frac{\partial f}{\partial t} + div(\rho \vect{u} f)=div\left(\displaystyle\frac{\alpha}{C}\, \grad f\right)+S
\end{equation}

\item[-] The coefficient $\alpha$ represents the sum of the
molecular and turbulent conductivities (depending on the models used), or
$\alpha=\alpha_m+\alpha_t$, with, for a turbulent viscosity based modeling,
$\displaystyle\alpha_t=C\,\frac{\mu_t}{\sigma_t}$, where $\sigma_t$ is the
turbulent Prandtl number\footnote{The turbulent Prandtl number is dimensionless and,
in some usual cases, taken as equal to $0.7$.}.

\item[-] The coefficient $C_p$ represents the specific heat, with unit $m^{2}\,s^{-2}\,K^{-1}=J\,kg^{-1}\,K^{ -1}$.

\item[-] We note $\lambda$ the thermal conductivity, with unit
$kg\,m\,s^{-3}\,K^{-1}=W\,m^{-1}\,K^{-1}$.

\item[-] It should be specified that $C=1$ for all variables except for the
temperature, for which $C=C_p$\footnote{More precisely, we have
$C=C_p$ for all {\it scalar variables} $f$ that we want to consider as the
temperature for the boundary conditions. These {\it scalar variables} can
be checked using the \var{is\_temperature} field keyword. By default this
value is set to 0 for all variables except for the temperature.}.
In the code, the value which the user must define is
$\displaystyle\frac{\alpha_m}{C}$ (if the property is constant, values are
initialized to \var{viscl0} for the velocity and \var{visls0} for {\it scalar variables};
If the property is variable, equivalent field values must be defined in the GUI or
in \fort{cs\_user\_physical\_properties.c}).

\item[-] For the variance of the fluctuations of a {\it scalar variable}, the
conductivity $\alpha$ and the coefficient $C$ are inherited from the
associated variable.

\end{itemize}


{\bf Simple Dirichlet condition}: for a simple Dirichlet condition, we
naturally obtain (special case of(\ref{Base_Condli_eq_fbord_condli})):
\begin{equation}
\underbrace{\ \ f_{b,int}\ \ }_{\text{boundary value used by the calculation}}
= \underbrace{\ \ f_{\text{\it real}}\ \ }_{\text{real value imposed on the boundary}}
\end{equation}
{\bf Other cases}: when the boundary condition is based on a flux,
this is of a diffusive flow \footnote{Indeed, the total outgoing flow of the domain is
given by the
sum of the convective flux (if the variable is actually convected)
and diffusive flux. However, for solid walls and symmetries, the mass flow is zero and
the condition reduced to a constraint on the diffusive flux. Moreover, for the
outlets (exiting mass flow), the boundary condition relates only to the
diffusive flux (often a homogeneous Neumann condition), convective flux
dependent on the upstream conditions (therefore it does not need
boundary conditions). Finally, at the inlets, a simple Dirichlet condition is most often
used, and the diffusive flux is deduced from it.}.
We then have:
\begin{equation}
\underbrace{\ \ \phi_{int}\ \ }_{\text{diffusive flux towards the internal domain}}
= \underbrace{\ \ \phi_{\text{\it real}}\ \ }_{\text{real diffusive flux imposed on the boundary}}
\end{equation}

The imposed real diffusive flux can be given
\begin{itemize}
\item [-] directly (Neumann condition), or
$\phi_{\text{\it r\'eel}}=\phi_{\text{\it imp,ext}}$ or
\item [-] implicitly deduced from two imposed values: an exterior value
$f_{imp,ext}$ and an exchange coefficient $h_{imp,ext}$
(generalized Dirichlet condition).
\end{itemize}

\vspace{1cm}
Depending on the type of condition (Dirichlet or Neumann) and assuming
the conservation of the flux normal to the boundary,
we can write (see figure \ref{Base_Condli_fig_flux_condli}):
\begin{equation}\label{Base_Condli_eq_flux_condli}
\begin{array}{l}
    \underbrace{h_{int}(f_{b,int}-f_{I'})}_{\phi_{int}}
  = \underbrace{h_{b}(f_{b,ext}-f_{I'})}_{\phi_{b}}
  = \left\{\begin{array}{ll}
    \underbrace{h_{imp,ext}(f_{imp,ext}-f_{b,ext})}_{\phi_{\text{\it imposed real}}} &\text{(condition of Dirichlet)}\\
    \underbrace{\phi_{\text{\it imp,ext}}}_{\phi_{\text{\it imposed real}}}
            &\text{(Neumann condition)}
           \end{array}\right.
\end{array}
\end{equation}

The ratio between the coefficient $h_{b}$ and the coefficient $h_{int}$ accounts for
the importance of crossing the near-boundary zone and is of particular importance in
the case of walls along which a boundary layer develops (whose properties are then taken
into account by $h_{b}$: see \fort{cs\_boundary\_conditions\_set\_coeffs\_turb}).
In the simpler framework considered here, we will limit ourselves to
case $h_{b}=h_{int}$ and $f_{b,ext}=f_{b,int}=f_{b}$.
The relation (\ref{Base_Condli_eq_flux_condli}) is then written:

\begin{equation}
\underbrace{h_{int}(f_{b}-f_{I'})}_{\phi_{int}}
  = \left\{\begin{array}{ll}
    \underbrace{h_{imp,ext}(f_{imp,ext}-f_{b})}_{\phi_{\text{\it imposed real}}} &\text{(condition of Dirichlet)}\\
    \underbrace{\phi_{\text{\it imp,ext}}}_{\phi_{\text{\it imposed real}}}
            &\text{(condition  Neumann)}
           \end{array}\right.
\end{equation}

By rearranging, we get the boundary value:
\begin{equation}\label{Base_Condli_eq_fbord_condli}
f_{b}
  = \left\{\begin{array}{cccccl}
    \displaystyle\frac{h_{imp,ext}}{h_{int}+h_{imp,ext}}&f_{imp,ext}&+&
    \displaystyle\frac{h_{int}}{h_{int}+h_{imp,ext}}    &f_{I'}
                         &\text{(Dirichlet condition)}\\
    \displaystyle\frac{1}{h_{int}}&\phi_{\text{\it imp,ext}}&+&
    \ &f_{I'}
            &\text{(Neumann condition)}
           \end{array}\right.
\end{equation}


{\bf Conclusion}: we will therefore note the boundary conditions
in the general form:
\begin{equation}
f_{b}=A_b + B_b\,f_{I'}
\end{equation}
with $A_b$ et $B_b$ defined according to the type of conditions:
\begin{equation}
\begin{array}{c}
\text{Dirichlet}\left\{\begin{array}{ll}
    A_b = &\displaystyle\frac{h_{imp,ext}}{h_{int}+h_{imp,ext}}f_{imp,ext}\\
    B_b = &\displaystyle\frac{h_{int}}{h_{int}+h_{imp,ext}}
                  \end{array}\right.
\text{\ \  Neumann}\left\{\begin{array}{ll}
    A_b = &\displaystyle\frac{1}{h_{int}}\phi_{\text{\it imp,ext}}\\
    B_b = &1
                  \end{array}\right.
\end{array}
\end{equation}

\newpage

{\bf Remarks }
\begin{itemize}
\item [-] The value $f_{I'}$ is calculated using the cell gradient of $f$,
that is: $f_{I'}=f_{I}+\vect{II'}\grad{f}_I$.
\item [-] The value of $h_{int}$ must yet be specified. This is a {\it numerical}, value,
having {\it a priori} no relation to a physical exchange coefficient,
which depends on the mode of calculation of the diffusive flux in the boundary cell.
Thus~$\displaystyle h_{int}=\displaystyle\frac{\alpha}{\overline{I'F}}$
(the unit is naturally deduced).
\item [-] We also repeat, because it can be a cause of error, that in the code, we have:
        \begin{itemize}
        \item [-] for the temperature $\alpha_m=\lambda$ et $C=C_p$
        \item [-] for the enthalpy      $\alpha_m=\displaystyle\frac{\lambda}{C_p}$ et $C=1$
        \end{itemize}
\end{itemize}

{\bf Examples of special cases}
\begin{itemize}
\item [-] In the case of a Dirichlet condition,
the user is therefore led to
provide two values: $f_{imp,ext}$ and $h_{imp,ext}$.
To obtain a simple Dirichlet condition (without exchange coefficient)
just impose $h_{imp,ext}=+\infty$. This is the most commonl case
(in practice, $h_{imp,ext}=cs\_math\_infinite\_r$ ).
\item [-] In the case of a Neumann condition, the user provides a single value
$\phi_{\text{\it imp,ext}}$ (zero for homogeneous Neumann conditions).

\end{itemize}

\etape{Output condition for pressure\vspace{0,3cm}}\label{Base_Condli_Sortie_Pression}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
We specify here condition applied to the
pressure for standard outlets.
It is necessary to impose a Dirichlet type condition (accompanied by a
homogeneous Neumann condition on the components of the velocity). We compute it
based on the values of the variable at the previous time step.

\begin{itemize}
\item [-] Reasoning on a simple configuration (a channel, with a flat output,
perpendicular to the flow), it can be assumed that the shape of the pressure
profiles taken on the surfaces parallel to the outlet are unchanged around
this outlet (hypothesis of an established flow, far from any
perturbation). In this situation, we can write
$\displaystyle\frac{\partial^2P}{\partial\vect{n}\partial\vect{\tau}_i}=0$
($\vect{n}$ is the vector normal to the outlet, $\vect{\tau}_i$ represents
a basis of the exit plane).

\item [-] If, moreover,
it can be assumed that the pressure gradient taken in the direction
perpendicular to the outlet faces is uniform in its neighborhood,
the profile to impose at the output (values $p_b$) can be deduced from the profile
taken on an upstream plane ($p_{upstream}$ values) by simply adding the constant
$R=d\,\grad{(p)}.\vect{n}$ (where $d$ is the distance between the upstream plane
and the outlet), or $p_b=p_{upstream}+R$ (the fact that $R$ is identical for
all outlet faces is important so we can eliminate it in the equation
(\ref{Base_Condli_eq_psortie_condli}) below).

\item [-] With the additional assumption that the points $I'$ relative to
the outlet faces are on a plane parallel to the output, we can use the
values at these points ($p_{I'}$) for upstream values, so
$p_{upstream}=p_{I'}=p_{I}+\vect{II'}.\grad{p}$.

\item [-] Furthermore, the
pressure being defined relative to a constant (incompressible flow)
one can fix its value arbitrarily at a point $A$
\footnote{by default, the center of mass of the outlet faces}) to $p_0$ (value fixed by
user, equal to \var{P0} and null by default),
and therefore shift the imposed profile at the output by adding:\\
$R_0=p_0-(p_{upstream,A}+R)=p_0-(p_{I',A}+R)$.

\item [-] So we finally obtain:
\begin{equation}\label{Base_Condli_eq_psortie_condli}
\begin{array}{lll}
p_b&=&p_{I'}+R+R_0\\
   &=&p_{I'}+R+p_0-(p_{I',A}+R)\\
   &=&p_{I'}+\underbrace{p_0-p_{I',A}}_{\text{constant value $R_1$}}\\
   &=&p_{I'}+R_1
\end{array}
\end{equation}
\end{itemize}
It is therefore noted that the pressure condition at the outlet is Dirichelet condition
whose values are equal to the pressure values (taken at previous time step) on
the plane upstream of the points $I'$ and readjusted to obtain \var{P0} in
an arbitrary exit point.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Remaining issues}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\etape{Representation of conditions by face value}
Although the method used allows simplicity and
homogeneity of treatment of all boundary conditions,
it is quite restrictive in the sense that a single value is not always sufficient to
represent the conditions to be applied when calculating different terms.

Thus, in $k-\varepsilon$ it was necessary, when calculating the boundary conditions of the wall,
to use two ($A_b$, $B_b$) pairs in order to take into account the conditions to apply*
for the calculation of the shear stress and those to be used when calculating the production
term (and a third set of coefficients would be necessary for
allow the treatment of the gradients intervening in the terms of gradient
transposed, in \fort{visecv}).

Perhaps it could be useful to set up a method
allowing to direclty use (at least in some strategic points of the code)
forces, stresses or fluxes, without requiring a general a face value computation.

\etape{Pressure outlet condition}
The outlet pressure condition translates to $p_f=p_{I'}+R1$ and the profile obtained
corresponds to the upstream profile taken at points $I'$ and readjusted to obtain $p_0$ at
an arbitrary point $A$. This type of condition is applied without precautions, but is not
always justified (a Dirichlet condition based on the calculated value directly at the boundary
faces might be more suitable).
The hypotheses are particularly faulty in the following cases:
\begin{itemize}
\item [-] the outlet is close to an area where the flow is not established
in space (or varies in time);
\item [-] the outlet is not perpendicular to the flow~;
\item [-] the pressure gradient in the direction normal to the outlet is not
the same for all output faces (in the case of multiple outputs, for example, the gradient
is probably not the same across all outputs);
\item [-] points $I'$ are not on a surface parallel to the output
(case of irregular meshes for example).
\end{itemize}

Moreover, in the absence of an outlet condition, it could perhaps
prove useful to set a reference value on a given cell
or to bring the average pressure back to a reference value (with the
shift of the spectrum, one ensures the invertibility of the matrix at each step
time, but we should also check whether the pressure is not likely to drift during
of the calculation).

\etape{Terms unaccounted for}
The current boundary conditions seem to cause issues when treating the transposed
velocity gradient term in the Navier-Stokes equations (handled as explicit term in the
time scheme). This term can be deactivated setting the \var{ivisse} keyword to $0$
in the velocity-pressure model.
